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FOREWORD 


The  research  effort  documented  in  this  report  was 
performed  under  the  Advanced  Technology  Torpedo  Thrust 
Program  and  funded  by  Office  of  Naval  Technology.  A  general 
purpose  program  has  been  developed  which  applies  the  concepts 
of  the  NSWC  Burn  Model  and  aenerates  a  tahulate^  data  set  . 
From  chis  raoulated  data  set,  any  one-,  two-,  or  three- 
dimensional  hydrocode  containing  a  reactive  mixture  scheme 
can  determine  a  new  burn  fraction  based  on  a  previous  burn 
fraction,  pressure,  and  time  step. 

The  illustrations  in  this  report  appear  after  the  body 
of  the  report. 

The  author  wishes  to  acknowledge  the  helpful  insights  of 
Dr.  Raafat  Guirguis  of  the  Naval  Surface  Warfare  Center  on 
how  the  NSWC  Burn  Model  works,  and  Drs .  Schittke  and  Feisler 
of  Indust rieanlagen-Betriebsgesellschaft  foi  their  suggestion 
to  develop  a  general  purpose  program  to  tabulate  the  data  and 
an  example  of  how  to  proceed  with  the  effort . 


Approved  by: 


WILLIAM  H.  BOHLI,  Acting 
Energetic  Materials  Division 
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SECTION  1 
INTRODUCTION 


This  report  describes  the  transformation  of  the  NSWC 
Burn  Model,  described  in  Reference  1,  into  a  general  purpose 
program  which  produces  a  simple  set  of  information.  This 
information  can  be  transformed  into  a  tabular  database  that 
can  be  used  by  any  one-,  two-,  or  three-dimensional  hydrocode 
which  contains  a  reactive  mixture  scheme.  The  term,  general 
purpose  program,  is  used  throughout  this  report  to  emphasis 
the  fact  that  the  computer  program  which  produces  the 
tabulated  data  set  is  separate  from  the  hydrocode  run  which 
uses  the  tabulated  data  set. 

The  basic  equations,  described  in  Reference  1,  which 
present  the  concepts  of  the  NSWC  Burn  Model  are  solved 
explicitly  by  a  general  purpose  program  in  a  Lagrangian 
coordinate  system.  These  equations  and  the  solution  method 
are  described  in  Section  2.  Relationships  which  prescribe 
the  necessary  parameters  for  a  hydrocode  computation  are 
shown  in  Section  3.  Section  4  presents  some  conclusions  and 
direction  of  future  work. 
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SECTION  2 

HOW  THE  NSWC  BURN  MODEL 
IS  USED  IN  A  GENERAL  PURPOSE  PROGRAM 


The  NSWC  Burn  Model  assumes  that  a  typical  pore  and 
binder  in  a  high  explosive  can  be  modeled  by  a  void  region 
inside  a  sphere.  Therefore,  it  is  referred  to  as  a  hollow 
sphere  model.  The  size  of  the  hollow  portion  of  the  sphere 
representing  the  void  and  binder  is  in  proportion  with  the 
global  percentages  of  the  component  explosive,  binder,  and 
pore.  The  solid  portion  of  the  hollow  sphere  located  between 
rif  the  minimum  radius  of  the  solid  portion  of  the  hollow 
sphere,  and  r0,  the  maximum  radius  of  the  solid  portion  of  the 
hollow  sphere,  is  composed  of  the  component  explosive. 

Figure  1-a  shows  a  typical  situation  in  a  real  component 
explosive.  Figure  1-b  shows  how  the  NSWC  Burn  Model  views  a 
typical  situation  in  a  real  component  explosive.  Figure  1-c 
shows  the  region  of  computation  of  the  NSWC  Burn  Model. 

The  NSWC  Burn  Model  divides  the  portion  of  the  sphere 
between  ri  and  r0  into  a  number  of  one-dimensional  subcells 
which  can  be  equally  or  geometrically  spaced.  In  general,  it 
seems  that  geometric  spacing  is  more  advantageous  since  the 
region  requiring  the  greatest  resolution  is  near  ri .  This  is 
due  to  the  fact  that  the  NSWC  Burn  Model  is  only  used  to 
determine  the  global  burn  fraction  until  a  certain  criteria 
is  reached,  indicating  that  surface  burning  has  become 
predominant  in  determining  the  burn  fraction.  Then  the 
global  burn  fraction  can  be  determined  according  to  a  surface 
burn  model. 
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It  is  assumed  in  the  general  purpose  program  presented 
here  that  the  temperature  throughout  the  hollow  portion  of 
the  sphere  representing  the  void  and  binder  is  constant.  The 
value  used  for  this  region  is  equal  to  the  value  of  the 
temperature  at  the  location  rx  initially,  and  at  the  location 
corresponding  to  the  inner  surface  of  the  solid  portion  of 
the  sphere  as  it  moves  into  the  solid  region  due  to  the 
conversion  of  solid  subcells  into  burned  gas  products  as  time 
progresses.  As  time  progresses,  the  small  portion  of  burned 
gas  products  over  all  the  subcells  is  collected  in  the  hollow 
portion  of  the  sphere  in  a  cumulative  manner. 

In  the  general  purpose  program  describing  the  NSWC  Burn 
Model,  each  subcell  is  bounded  by  two  location  values,  the 
difference  of  which  is  the  subcell  width.  The  local  burn 
fraction  and  temperature  are  found  at  the  center  of  each 
subcell.  The  velocity  is  found  at  the  two  locations  which 
bound  the  subcell  region.  Initially,  it  is  assumed  that  the 
local  burn  fraction  is  zero,  the  temperature  is  the  ambient 
temperature,  and  the  velocity  is  zero  for  all  subcells. 

Initially,  a  pressure  due  to  the  surroundings,  pD,  is 
applied  at  r0  and  the  pressure  of  the  gas  in  the  hollow 
portion,  pg,  of  the  sphere  is  ambient  pressure  or  some 
specified  condition  and  applied  at  r^  As  the  pressure,  pQ, 
is  applied  to  the  outer  surface  of  the  hollow  sphere  model, 
it  is  assumed  that  there  is  a  global  mechanical  deformation 
in  progress  near  the  pore,  and  that  the  material  behavior  is 
such  that  the  stress  is  dependent  on  strain  and  strain  rate. 
This  results  in  an  increase  in  pg .  The  processes  which  occur 

are  described  in  greater  detail  in  Reference  1 .  Two 
assumptions  are  made,  the  elastic  portion  of  mechanical 
deformation  can  be  ignored  and  the  bulk  compressibility  of 
the  component  explosive  can  be  neglected.  This  greatly 
simplifies  the  equations  needed  to  describe  the  NSWC  Burn 
Model  for  a  hollow  sphere. 
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Also,  the  location  of  the  subcells  in  the  hollow  sphere 
model  moves  toward  the  center  of  the  sphere  as  the  pressure, 
pD,  is  applied  to  the  outer  surface  of  the  hollow  sphere 
model.  Likewise,  as  the  component  explosive  begins  to  burn, 
the  gas  pressure,  p  ,  in  the  hollow  portion  of  the  sphere 
begins  to  increase,  resulting  in  a  force  being  applied  to  the 
surface  located  at  rA  and  directed  toward  the  outer  surface  of 
the  sphere.  This  results  in  additional  changes  to  subcell 
boundary  locations.  The  equation  used  to  determine  the 
velocity,  v(r,x,t),  of  the  boundaries  associated  with  each 
subcell  is  the  following, 1 


7  (Po  ~  Pg  ~  2 V~3  k  ln(r0/ri)) 
2  (r0~3  -  r^3)  r2  k 


(1) 


The  variables  used  in  Equation  {1)  are  described  below. 


The  basic  equations  used  in  the  general  purpose  program 
to  describe  a  hollow  sphere  model  representative  of  the  NSWC 
Burn  Model  in  which  an  Arrhenius  reaction  is  applied  to  the 
hot  spot  temperature  which  is  modified  by  thermal  diffusion, 
are  the  following,1 


dX (x, t ) 
dt 


err 


dA ( r , x , t ) 

_d£ _ 


4  71  r2 


~  n  (rQ3  -  rq3) 


dr  , 


dA ( r , x ,  t ) 
dt 


:i  -  A(r,x,t))  Z  exp - 


T  (r ,  x, 1 1 


(2) 


(3) 


T (r, x, t) 


To (r, x, t) 


+ 


^  ( r ,  x,  t )  dt 
dt 


and 


(4) 
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dT  (r,  x,  t )  =  2.25  y  (po  -  pg  -  2  j~3  k  ln(r0/rj))2 
dt  p  Cp  (r^3  -  r0-^|  k 

+  1  1  d  lr2  dT  (r, x, t)  |  +  _Q_  dA  (r,  x,  t )  (5) 

p  Cp  r2  dr  \  dr  I  Cp  dt 

where  r  is  the  location  of  the  subcell  being  evaluated  at 

time,  t;  X(x,t)  is  the  global  burn  fraction  of  the  system; 

A(r,x,t)  is  the  local  burn  fraction  of  an  individual  suocell 

in  the  hollow  sphere  model;  and  T(r,x,t)  is  the  temperature 

of  an  individual  subcell  in  the  hollow  sphere  model. 

Equations  (3)  through  (5)  are  solved  explicitly  in  the 

Lagrangian  subcell  coordinate  system.  The  variable  x 

indicates  the  global  computational  cell  for  which  the  hollow 

sphere  model  is  being  computed.  For  example,  in  a  three- 

dimensional  Cartesian  coordinate  system,  x  would  be  a 

function  describing  a  unique  X,  Y,  and  Z  location  for  the 

global  computational  cell.  In  the  general  purpose  program 

described  in  this  paper,  the  hollow  sphere  model  is  evaluated 

for  one  computational  cell.  Therefore,  x  is  set  to  1. 

T0(r,x,t)  is  the  initial  temperature  of  an  individual  subcell 

in  the  hollow  sphere  model,  Z  is  the  pre-exponential  factor 

in  the  Arrhenius  kinetics,  T  is  the  activation  temperature,  p 

* 

is  the  density,  k  is  the  thermal  conductivity,  Cp  is  the  heat 
capacity,  Q  is  the  heat  of  reaction  of  the  explosive,  and  y  is 
a  number  which  affects  how  fast  the  global  burn  fraction 
ri^es  but  not  the  asymptote  to  which  the  global  burn  fraction 
approaches.  The  values  for  the  constants  used  to  model 
PBX-9404*  are  shown  in  Table  1  on  the  next  page. 

Due  to  the  nature  of  Equation  (4)  for  dT (r, x, t ) /dt ,  it 
is  necessary  to  determine  the  quantity, 

-r_  -1-  -1  (r*  k-  aTlr-x-t)|  ,  (6) 

p  Cp  r2  dr  \  dr  I 
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TABLE  1.  PBX-9404  PARAMETERS  USED  IN  THE  NSWC  BURN  MODEL 


z 

5.0  x  1013 

(lsec 

T* 

26500 

°K 

T0 

300 

°K 

P 

1 . 84 

,  3 

g/cm 

CP 

L  '• 

1 

o 

f-H 

X 

t — i 

2  2  c 

cm  /|lsec  /  K 

k 

8.0  x  10~5 

Mbar 

Q 

5.439  x  10~2 

cm2/|lsec" 

k* 

8  x  10"14 

cm/|isec/°K 

ri 

0.0039 

cm 

ro 

0.01 

cm 

for  both  boundaries  of  a  subcell.  This  takes  into  account 
both  the  temperature  transferred  from  the  neighboring  subcell 
on  one  side  of  a  particular  subcell  and  the  temperature 
transferred  to  the  neighboring  subcell  on  the  other  side  of  a 
particular  subcell. 

The  general  purpose  program  which  uses  Equations  (1) 
through  (6)  is  shown  in  Appendix  A.  When  running  this 
program,  the  user  is  required  to  enter  the  time  step  for 
which  the  equations  will  be  sol /ed  explicitly  and  the  initial 
pressure,  p0,  applied  to  the  hollow  sphere.  The  user  may 
spec.*  fy,  in  the  form  of  data  statements  in  the  general 
purpose  program,  the  number  of  subcells  in  the  component 
explosive  portion  of  the  sphere  and,  if  geometric  cell  size 
expansion  is  used,  the  initial  cell  size  and  expansion  factor 
from  one  cell  to  the  next.  The  parameters  shown  in  Table  1 
are  included  in  the  form  of  data  statements  in  the  general 
purpose  program. 
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The  general  purpose  program  iterates  twice  through  the 
section  which  dete~mines  dX(x,t)/dt,  dA (r, x, t) /dt,  T(r,x,t), 
and  dT(r,x,t)/dt  to  converge  to  the  solution.  Iterations 
greater  than  two  do  not  add  significant  improvements  in 
convergence.  The  pressure  of  the  gj  s  is  determined  by  usin._, 
a  JWL  ( Jones-Wilkens-Lee)  Equation  of  State  to  model  the 
component  explosive  gt.s  prouucts.  The  JWL  Equation  of  State 
has  the  following  form, 

p  =  A  (1 - — — )  exp  (-R;  V'  +  B  (1  -  — )exp(-R2  V)  +  ^  E  (7) 

Rj  V  R2  V  V 

where  V,  the  relative  volume,  equals  v/v0  or  po/p.  The  values 
for  the  JWL  Equation  of  State  parameters  used  to  model 
PBX-94042  are  shown  in  Ta»_le  2. 

TABLE  2.  JWL  EQUATION  OF  STATE  PARAMETERS  USED  ~OR  PBX-9404 


p 

1.84 

,  3 

g/cm 

Pcj 

0.37 

Mbar 

D 

0.88 

cm/psec 

Ec 

0.102 

Mbar-cc/cc 

Ri 

4  .  6 

R? 

1.3 

(h 

0.38 

A 

8 . 524 

Mbar 

B 

0 . 1802 

Mbar 

C 

0.0121 

Mbar 

The  relative  volume,  V,  is  found  by  the  fo1 lowing 
method.  First  one  multiplies  the  volume  between  r-  and  r.  by 
the  globax  burn  fraction.  This  results  in  a  quantity,  v.,, 
which  takes  into  account  the  relative  volume  of  the  high 
explosive  gas  products  from  the  component  explosive.  If  v,  is 
less  than  l.E-11,  it  is  set  equal  to  l.E-11.  The  current 
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volume  of  the  hollaw  sphere,  /hs,  and  the  volume  of  the 
binder,  vD,  are  also  included  in  the  determinat’on  of  the 
relative  volume  which  is  given  by  the  following  expression, 


V 


vhs  +  vq  -  vb 


v 


g 


(8) 


Initially,  it  was  desirable  to  create  a  model  with 
constants  that  were  solely  determined  by  the  physical 
parameters  of  the  component  explosive.  However,  the  use  of 
the  JWL  Equation  of  State  to  model  the  component  explosive 
gas  products  introduces  another  set  of  unknowns  which  are 
usually  determined  through  experiments  such  as  cylinder  tests 
or  similitude  data.  In  future  applications  of  the  general 
purpose  program  for  high  explosives  that  do  not  have  JWL 
Equation  of  State  parameters,  it  will  be  assumed  that  a 
generic  set  of  JWL  Equation  o*  State  parameters  will  be 
sufficient  to  model  the  gas  products. 
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SECTION  3 

HOW  TO  APPLY  THE  GENERAL  PURPOSE 
PROGRAM  RESULTS  TO  A  HYDRODYNAMIC  COMPUTATION 


The  general  purpose  program  developed  in  the  previous 
section  for  the  NSWC  Burn  Model  was  designed  to  evaluate  one 
hot  spot  formation  caused  by  a  constant  pressure  applied  at 
the  outer  surface  of  the  hollow  sphere  for  the  duration  of 
the  run.  In  a  one-,  two-,  or  three-dimensional  hydrodynamic 
computation,  there  are  many  cells;  each  cell  experiences 
varying  pressure  levels  throughout  the  computation.  In  this 
section,  a  scheme  for  using  the  results  of  the  general 
purpose  program  describing  the  NSWC  Burn  Model  will  be 
presented . 

The  first  characteristic  which  needs  to  be  defined  is 
the  transition  from  exponential  growth  of  the  reaction  due  to 
mechanical  deformation  to  a  slower  growth  of  the  reaction  due 
to  surface  burning.  This  exponential  growth  of  the  reaction 
due  to  mechanical  deformation  can  be  thought  of  as  the 
ignition  phase  of  a  hot  spot.  The  equations  and  general 
purpose  program  presented  in  Section  2  describe  this 
phenomena.  The  growth  phase  in  which  a  transition  to  surface 
burning  occurs  is  described  with  a  different,  simpler 
expression  for  the  global  burn  fraction,  such  as  a  surface 
burn  model.  This  expression  is  described  in  more  detail  in 
Reference  1  . 

After  evaluating  how  different  system  parameters  changed 
as  the  growth  of  the  reaction  slowed,  it  was  observed  that 
the  quantity, 
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A.(x,t)  *  ^T(j,x,t)  *  dr  (  j  ,  x,  t ) 

_ j^i _ 

n  *  (rG  -  ri)  p0 

changed  sign,  from  a  negative  to  a  positive  value,  when  the 
exponential  growth  of  the  reaction  due  to  mechanical 
deformation  slowed  for  every  applied  pressure  (from  5  to  200 
kilobar) .  The  quantity  represented  by  Equation  (9)  is 
independent  of  the  time  step  step  selected.  In  Equation  (9), 
X(x,t)  is  the  global  burn  fraction,  T(j,x,t)  is  the 
temperature  in  the  subcell,  n  is  the  number  of  subcells 
describing  the  component  explosive,  dr(j,x,t)  is  the  width  of 
the  subcell,  r0  and  ^  are  the  current  values  of  the  solid 
portion  of  the  sphere  which  bound  the  component  explosive 
region  in  the  model,  and  pQ  is  the  applied  pressure  as 
described  in  the  previous  section.  When  the  quantity 
represented  by  Equation  (9)  changed  sign  in  the  general 
purpose  program,  it  was  assumed  that  the  global  burn  fraction 
reached  the  transition  point,  Xhs,  at  time,  ths,  where  surface 
burning  becomes  the  predominant  mechanism. 

When  running  a  one-,  two-,  or  three-dimensional 
hydrodynamic  computation,  another  criteria,  based  on  the 
computation,  is  used  to  determine  when  the  transition  from 
ignition  to  growth  should  be  made.  This  transition  to  growth 
is  based  upon  a  minimum  burn  rate  computed  from  both  the 
ignition  model  and  the  growth  model.  This  method  of 
determining  transition  (presented  above)  provides  enough 
information  to  reach  and  go  slightly  beyond  the  criteria  of  a 
hydrodynamic  computation. 

Important  consideration  also  needs  to  be  given  when 
determining  the  minimum  pressure  below  which  a  burn  fraction 
will  not  accumulate  in  hydrodynamic  computations.  The 
general  purpose  program  was  run  starting  at  a  pressure  of 
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200  kilobar.  In  each  successive  run,  the  pressure  was 
reduced  by  1  kilobar.  It  was  determined  that  5  kilobars  is 
the  lowest  pressure  that  produced  a  smooth  transition  of  X„s 
from  one  pressure  to  the  next.  At  6  kilobars,  the  value  of 
Xhs  was  4.92e-3;  at  5  kiloba.LS,  the  value  of  was  4.15e-3; 
and  at  4  kilobars,  the  value  of  Xhs  was  1.22e-8.  Therefore, 
the  lowest  acceptable  pressure  has  been  selected  as 
5  kilobars.  The  data  generated  at  1  kilobar  intervals 
includes  the  global  burn  fraction  at  transition,  Xrs,  and  the 
associated  time  at  which  this  occurs,  ths  for  each  pressure. 

For  pressures  of  5  kilobars,  pressures  of  10  to 
100  kilobars  at  10-kilobar  intervals,  and  pressures  of  100  to 
200  kilobars  at  20-kilobar  intervals,  the  burn  fraction 
versus  time  data  is  written  out  for  every  time  step.  This 
produced  a  large  amount  of  data  which  was  reduced  to  be 
easily  handled  by  a  hydrodynamic  computation.  The  program 
which  reduces  the  data,  SLIM,  is  shown  in  Appendix  B.  In 
addition  to  data  reduction,  SLIM  also  normalizes  the  data, 
determines  the  length  of  the  line  resulting  from  the 
normalized  data  points,  and  determines  the  coordinates  of 
points  which  divide  the  normalized  line  up  into  400  equally 
spaced  segments  based  on  the  original  normalized  data  points. 
This  allows  the  hydrodynamic  computation  to  interpolate 
between  two  ranges  of  pressure  to  generate  a  normalized  burn 
fraction  versus  time  curve  for  any  pressure.  It  also  can  be 
used  to  find  the  new  burn  fraction. 

After  much  investigation,  the  following  normalization 
scheme  was  found, 


where  Y  is  0.15  when  X  is  the  time,  Y  is  0.05  when  X  is  the 
burn  fraction,  P  is  the  pressure  in  the  cell,  and  Pc-  is 
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the  Chapman-Jouguet  pressure  of  the  high  explosive.  This 
normalization  scheme  seemed  to  work  fairly  well  in  creating 
curves  from  which  one  could  interpolate.  The  normalized 
curves  for  5  kilobars,  10  through  100  kilobars  at  10-kilobar 
intervals,  and  100  through  200  kilobars  at  20-kilobai 
intervals  are  shown  in  Figure  2. 

In  order  to  use  this  information  in  a  hydrodynamic 
computation,  one  tabulates  the  normalized  curves  in  the  form 
of  data  statements  so  that  a  hydrodynamic  computation  can 
determine  the  required  information.  The  two  programs,  MAKEHS 
and  MAKELAM,  used  to  create  the  necessary  data  statements  are 
shown  in  Appendix  C.  A  sample  subroutine  which  uses  this 
information  is  shown  in  Appendix  D.  Currently,  this 
subroutine  is  used  in  the  Eulerian  hydrocode  DYSMAS . 

The  first  time  a  cell  experiences  a  pressure  greater 
than  a  predetermined  minimum,  such  as  5  kilobars  for  the  case 
of  PBX-9404,  the  value  of  the  normalized  time  is  determined 
using  the  time  step  as  t,  the  current  pressure  in  the  cell  as 
pQ,  ths  based  on  current  p0,  and  the  Chapman-Jouguet  pressure 
of  the  component  explosive.  From  this,  the  associated  burn 
fraction  and  burn  rate  can  be  determined.  On  subsequent 
cycles,  the  previous  burn  fraction  is  used  to  determine  an 
old  normalized  time  based  on  the  current  value  of  p0 .  Then 
the  current  time  step  is  used  to  determine  a  new  normalized 
time  and  associated  burn  fraction  ana  burn  iate.  Using  this 
method,  it  is  quite  feasible  to  model  the  ignition  phase  with 
the  NSWC  Burn  Model  in  conjunction  with  a  simple  surface  burn 
growth  phase  model,  and  to  run  one-,  two-,  or  three- 
dimensional  hydrodynamic  computations. 
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SECTION  4 

CONCLUSIONS  AND  FUTURE  WORK 


The  development  of  a  general  purpose  program  to 
incorporate  the  concepts  of  the  NSWC  Burn  Model  into  a 
tabular  data  set  have  been  successful.  Future  work  will 
involve  comparing  results  from  inclusion  of  the  tabulated 
data  set  into  various  hydrococies  to  results  shown  in 
Reference  1.  To  complete  this  work,  constants  in  the  surface 
burn  model  may  need  to  be  modified  to  match  the  results  shown 
in  Reference  1. 

Future  work  also  will  involve  (1)  verifying  the 
assumption  (discussed  in  Section  2)  that  a  generic  JWL 
Equation  of  State  is  sufficient  because  the  parameters  of  the 
NSWC  Burn  Model  will  compensate  as  needed,  and  (2) 
investigating  whether  the  criteria  for  transition  from 
ignition  to  growth  phase  will  always  be  the  quantity,  shown 
in  Equation  (9),  changing  sign  for  different  explosive 
components  and/or  similar  explosive  components  with  varying 
particle  sizes,  initial  porosities,  initial  temperatures, 
and/or  initial  mass  distributions  between  the  component 
explosive  and  binder. 

Planned  improvements  to  the  general  purpose  program 
include  (1)  implementing  a  numerical  expression  for  the  bulk 
compressibility  of  the  binder  to  take  into  account  the  energy 
loss  of  compressing  the  binder,  and  (2)  investigating  a 
method  of  including  shear-induced  hot  spot  temperature 
effects . 
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COMPONENT  EXPLOSIVE 

BINDER 

PORE  (VOID) 

AREAS  OF  HIGH  MECHANICAL 
DEFORMATION  (HOT  SPOTS) 


(a)  real  situation  in  composite  explosive 


i=l  to  n 


[Reg  ion  of  I 
Computation 


(b)  one-dimensional  approach 
of  the  NSWC  Burn  Model 


(c)  schematic  of  one-dimensional 
computational  region  of  the 
NSWC  Burn  Model 


FIGURE  1. 


THE  NSWC  BURN  MODEL 


200  Kbar 
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APPENDIX  A 

GENERAL  PURPOSE  PROGRAM  WHICH 
IMPLEMENTS  THE  CONCEPTS  OF  THE  NSWC  BURN  MODEL 


c  3/4/91  version 
program  test 

implicit  double  precision (a-h, o-z) 

dimension  temp ( 450 ) , clam (450) ,drrl3 (450) , rm(450) , 

1  tburn (450) , dslmax ( 450 )  ,  tempold ( 450 ) , rmold ( 450 ) 
dimension  rbd (451) , rbd3 (451) , dr (450) 
data  m/ 415/, fact in /l. 01/, cell in /0. 000001/ 
data  z/5.el3/,  tstar/2 . 65e4/,  tambi/300./,  rho/1.84/ 
data  cp/1.4e-5/,  bulk/8. e-5/,  q/5.439e-2/,  stark/ 8 . e-1 4 / 
data  ri/0.0039/,  ro/0.01/,  thrd/0 . 33333333/ 
data  ajwl/8. 524011/,  bjwl/0 . 1801812/,  omjwi/U .38/ 
data  rljwl/4.6/,  r2jwl/1.3/,  ojwl/0 . 012066373/ 
write  (6,*)  'input  time  step' 
read (5, * )  dto 
compold=10000 . 
totlen=ro-ri 
gam=0 . 025 
po=25 .e-3 

type  *, 'Enter  initial  pressure' 

read(*, *) po 

ms=l+ (m-1) /10 

mpl=m+l 

nc=0 

tt=0 . 

c  hot  spots  will  not  form  if  it  takes  longer  than  1  microseconds 
tl  =  l. 
dt=dto 
ib=l 
istr=l 
rg=ri 

rg3=rg*rg*rg 

rbd(l)=ri 

rbd3 ( 1 ) =rbd ( 1 ) *  *3 

ro3=ro*  *3 

ri3=rbd3 (1) 

ro3mri3=ro3-ri3 

fako=2 . 25*gam/ ( rho*cp*bulk* (1 . / ri3-l . /ro3) **2) 

fakr=3 . / (ro3mri3) 

fakp=po-3 ,464102*bulk*dlog(ro/ri) 

f akq=q/ cp 

f aks=l . /cp/ rho 

fakv=gam/2 . /bulk 

prev  =  cellin/f actin 
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sum  =  ri 
do  i=l/m 
ipl=i+l 

prev=f act in *p rev 
dr  (i) =prev 

if (sum+prev. It . ro) sum=sum+prev 
if (i .eq.m) dr (i) =ro-sum 
temp (i) =tambi 
tempold(i) =tambi 
clam(i) =1 .e-30 
tburn ( i) =0  . 
dslmax (i) =0 . 
rbd(ipl) =dr (i) +rbd(i) 
rbd3 ( ipl ) =rbd ( ipl ) **3 
drrl3 (i) =rbd3 (ipl) -rbd3 (i) 
rm (i) =rbd (ipl) -0 . 5*dr (i) 
rmold(i) =rm(i) 
end  do 
10  continue 

determine  volume  of  gas  if  not  first  time  through 

vo=slam* (rbd3 (mpl) -rbd3 (1) ) 
if (vo . It . 1 .d-11) vo=l .d-11 
vbold=vb 

vb=0 . 01*rbd3 (mpl) 
spvolg= ( rbd3 (1) +vo-vb) /vo 
pgold=pg 

pg=a jwl*dexp (-rl jwl*spvolg) +b jwl*dexp (-r2 jwl*spvolg> + 
1  c jwl*spvolg** (-omjwl-1) 
if (istr ,eq.l)pgbaseline=pg 
if (istr .eq. 1) pg=l . 01323e-6 
if (istr . ne . 1) pg=pg+l . 01323e-6-pgbaseline 

determine  new  radial  positions  based  on  velocity 

if (istr.ne.l) then 
riold=rbd(l) 
roold=rbd (mpl) 
do  i=l,mpl 

v=f akv* (fakp-pg) / (l/ro3-l/ri3) /rbd(i) **2 
rbd (i) =rbd (i) +v*dt 
rbd3 (i)=rbd(i) **3 
end  do 
ro=rbd (mpl ) 
ri=rbd(l) 
totlen=ro-ri 
ro3=rbd3 (mpl ) 
ri3=rbd3 ( 1 ) 
ro3mri3=ro3-ri3 
do  i=l,m 

dr (i) =rbd (i+1) -rbd (i) 
rm ( i ) =rbd ( i ) +dr ( i ) / 2 . 
end  do 

f akp=po-3 . 464l02*bulk*dlog(ro/ri) 
f ako=2 . 25*qam/ ( rho*cp*bulk* (1. /ri3-l . /ro3) **2) 
f akr=3 . / (ro3mri3) 


A- 2 


NAVSWC  TR  90-364 


endif 
tt=tt+dt 
dslam=0 . 0 
tempave=0 . 0 
dtempave=0 . 0 
do  i=l,m 

if (clam(i) .eq. 1 . dO) goto  55 
r=rm (i) 

r2dr=r*r*dr (i) 

clamold=clam(i) 

dslamold=dslam 

if ( i . ne . 1 . and . i . ne . m) then 

deltempr= (tempo Id (i+1) -tempo Id (i) ) / (rmold (i+1 ) -rmold ( i) ) 
deltempl= (tempo Id (i) -tempo Id ( i— 1 ) ) / (rmold (i) -rmold ( i— 1 ) ) 
elseif (i .eq. 1) then 

deltempr= (tempold (i+1 ) -tempold (i) ) / (rmold (i+1) -rmold (i) ) 
deltempl=0 . 0 
elseif (i .eq.m) then 
deltempr=0 . 0 

deltempl= (tempold (i) -tempold (i-1) ) / (rmold (i) -rmold (i-1 ) ) 
endif 

subdtempr=deltempr* ( (rm(i+l) +rm(i) ) /2 . ) **2 
subdtempl=deltempl* ( (rm(i) +rm(i-l) ) / 2 . ) **2 
subdtemp=f aks*stark* (subdtempr-subdtempl) /r2dr 
do  54  jj=l,2 

tdt=t star /temp (i) 

if (tdt . gt . 100 . ) type  *, ' t star /temp (i) 1 , tdt 
if (tdt .gt . 100 . ) tdt=100 . 

dclam= (1 . -clam(i) ) *z*dexp (-tstar/temp (i) ) 
dtemp=fako* (fakp-pg) **2/r**6+f akq*dclam+subdtemp 
clam(i) =dminl (1 .dO, clamold+dclam*dt ) 
if (clam(i) . eq . 1 . dO ) then 
dclam= (1 . -clamold) /dt 

dtemp=fako* (fakp-pg) **2/ r**6+f akq*dclam+subdtempr+subtempl 
endif 

temp (i) =tempold (i) +dtemp*dt 
dslam=dslamold+f akr*dclam*r2dr 
if (clam(i) .eq. 1 ,d0) goto  55 

54  continue 

55  continue 
tempave=tempa ve+temp ( i ) *dr ( i ) 
dtempave=dtempave+dtemp*dr (i) 
end  do 

istr=0 

slam=slam+dslam*dt 
do  56  i=l,m 

tempold (i) =temp ( i) 
rmold (i) =rm(i) 

56  continue 

rg3=rbd3 (ib) +clam(ib) *drrl3 (ib) 
rg=rg3**thrd 

if (clam(ib) .ge . 0 . 999999d0)  then 
tburn (ib) =tt 
dslmax (ib) =dslam 
ib=ib+l 
end  if 

nc=mod (nc+1, 10) 
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mss=5 

mdum=25 

comp=tempave/m/tot len*slam/po 
if (comp .eq . 0 . ) comp=compold-l . 
comp=dlogl0 (comp) 

c  hot  spot  lambda  is  reached  according  to  criteria  on  next  line 

if (comp . gt . 0 . 0 ) goto  35 

compold=comp 

slamold=slam 

if (po . eq. 0 . 0 1 . or . po . eq. 0 . 02 . or . po . eq ,0.03.or.po.eq.0.04 

1  .or.po.eq.0.05.or.po.eq.0.06.or.po.eq.0.07 

2  . or . po . eq .0.08.or.po.eq.0.09.or.po. eq .0.10 

3  . or . po . eq . 0 . 1 1 . or . po . eq .0.12.or.po.eq.0.13 

4  . or . po . eq. 0.14.or.po.eq.0.15.or.po.eq.0.005) 

1  write(8, ' (lx, Ip, 5el3 . 5) ') slam, dslam, comp, tt , pg 
if (tt . It . tl . and . clam (m)  .  It .  1 .  )  goto  10 
c  the  next  line  is  designed  to  crash  the  batch  job  when  t  hot  spot 
greater  than 
c  1  microsecond. 

if(tt.ge.tl)acrash=l./0. 

35  continue 

open (22,file='hs.info',status=' unknown  ' ) 
do  838  ids-1, 10000 

read (22, *, end=839, err=839)a,b,c 

838  continue 

839  continue 

write (22 ,  ’  (7 (lx,el5.?) )  ’  ) 

1  po, tt-dt, slamold, pgold, riold, roold, vbold 

close  (22) 

type  *, 't  hot  spot  =  ', tt-dt,'  slam  =  ', slamold 

stop 

end 
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APPENDIX  B 

PROGRAM  TO  REDUCE  AND  NORMALIZE 
OUTPUT  FROM  THE  GENERAL  PURPOSE  PROGRAM 


c  compile  with  for/g_float 

implicit  double  precision (a-h, o-z ) 
dimension  asave(lOOO) , dsave (1000) 
character*14  afile 
character*10  bfile 
aold=1000000 . 

:ype  *, 'enter  file  name  to  read1 
read ( * , 100 ) afile 
100  format (al 4 ) 

bfile (1 : l)=afile  <1 : 1) 
bfile(2:10)=a file (6:14) 
open (82,file=afile,status='old‘) 
open (84, f ile=bfile, status= ’new' ) 
inum=0 

do  120  i=l,  1  0000C 0 

read (82 , * ,  end=999) a, b, c, d, e 

anew=dlogl0 (a) 

if (dabs (anew-aold) . gt . . 1 . or . 

1  (a. gt. 0.0001. and. dabs (anew-aold) .gt . . 01 ) ) then 

inum=inum+l 
if ( inum. gt . 1000 ) stop 
asave ( inum) =a 
asave (inum) =d 
ipr~l 
aold=anew 
else 
ipr=0 
endif 

120  continue 
999  continue 

if ( ipr . eq . 0 ) then 
inum=i num+1 
if ( inum. gt . 1000) stop 
asave ( inum) =a 
dsa ve ( inum) =d 
endif 

type  *,'what  is  pressure, pc j ' 
read ( * , * ) pres , pc  j 
quanx= (pc j/pres) *  * . 15 
quany= (pc  j /pres ) *  * .05 
c  convert  to  scaled  quantities 
do  200  iic=l,inum 

asave ( iic) =dlogl0 ( (asave ( iic) /asave ( inum) ) *  *quany) 
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dsave (iic) =dlogl0 ( (dsave (iic) /dsave (inum) ) **quanx) 

200  continue 

c  find  total  length 
inumml=inum-l 
sumlen=0 . 

do  210  iic=l , inurranl 

sumlen=sumlen+dsqrt ( (asave(iic+l)-asave(iic) ) **2+ 

1  (dsave (iic+1) -dsave (iic) ) **2) 

210  continue 

icar=2 

icarml=icar-l 
dellen=sumlen/400 . 
curlen=0 . 0 
ruf len=0 . 0 
if lr=l 

write (94, 431) dsave (1) ,asave(l) 
do  220  i=l, 399 

curlen=curlen+dellen 
230  if ( if lr . eq . 1 ) ruf len=ruf lent 

1  dsqrt ( (asave (icar) -asave (icarml) ) **2+ 

2  (dsave (icar) -dsave (icarml) ) **2) 
if (curlen .gt .  ruf 1 en) then 

icarml=icar 
icar=icar+l 
if lr=l 
goto  230 
endif 
if lr=0 

dbtp=dsqrt ( (a  lave (icar) -asave (icarml )  )  **2  + 

2  (dsave (icar) -dsave (icarml ) ) *  *  2 ) 

ari=dbtp- (ruf len-curlen) 

aval=ari* (asave (icar) -asave (icarml) ) /dbtp+asave (icarml) 
dval=ari* (dsave (icar) -dsave (icarml) ) /dbtp+dsave (icarml) 
write (84, 431) dval, aval 
220  continue 

431  format (lx, fll .7, 4x, fll.7) 

write(84,431)dsave( inum) , asave ( inum) 

stop 

end 
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APPENDIX  C 

PROGRAMS  TO  CREATE  FORTRAN  DATA 
STATEMENTS  FROM  THE  TABULATED  DATA  SET 
FOR  INCLUSION  INTO  A  SUBROUTINE  CALLED  BY 
ANY  ONE-,  TWO-,  OR  THREE-DIMENSIONAL  HYDROCODE 


program  makehs 
parameter  iam=500 

dimension  po (iam) , ths (iam) , alamhs (iam) ,pg (iam) , ri (iam) , 

1  ro (iam) , vb (iam) 

open (81,file='hs.info,,status='old') 
icou=0 

do  100  i=l, 1000 

read ( 81 , * , end=999) po (i) , ths (i) , alamhs ii) , pg (i) ,  ri (i)  , 

1  ro  (i)  ,  vb  (i) 

icou=icou+l 
100  continue 

999  continue 

open (82 , f ile= 'hsdata . f or ' , status= ' new ' ) 
write (82, 110) icou, icou, icou, icou, icou 
110  format (6x, 'dimension  po ( ' , i3, ' ) , ri ( ' , i3, ’ ) , ro ( ' , i3, 

1  ’),',/,5x,'l  alamhs (', i3, '), ths (’, i3, ')' ) 

write (82, 115) icou 

115  format (6x,  'data  nrs/ ' , i3,  '  /  ' ) 

do  200  i=l, icou 

write(82,120)i,ri(i) , i, ro (i) , i, po (i) 

120  format ( 6x, ' data  ri  { ' , i3, ' ) / '  ,el3 .7, 

1  ' /, ro  ( ’ , i3,  ' ) / ' , el3 . 7, '/,po(',i3, ' ) / ' , f 5 . 3, '/') 

200  continue 

do  210  i=l, icou 

write(82,13rUi.  alamhs  (i)  ,  i,  ths  (i) 

130  format ( 6x,  ' data  alamhs (', i3,  ')/', el3 . 7,  ' /,  ths ( ' , i3,  ' ) / ' , 
)  el3 .7,  '  /  '  ) 

210  continue 

stop 
end 
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program  makelam 

implicit  double  precision  (a-h,o-z) 

character*10  filename 

parameter  iam=1000 

dimension  adum ( iam) , tdum ( iam) 

numf il=l 6 

pc j= . 370 

f i lename ( 1 : 2 )  =  ' f  8 ' 
filename (7:10)=' . dat 1 

open ( 82 , f ile= ' lamf ile .for',status=' new ' ) 

do  90  iiu=l,numfil 

if ( iiu . eq . 12 ) isep=iiu 

if ( iiu . gt . 12 ) isep=isep+2 

if (iiu.eq.16) filename (3 : 6) = 'p200 ' 

if (iiu .eq. 15) filename (3:6)='pl80' 

if (iiu . eq . 14 ) filename (3:6)='pl60' 

if (iiu . eq . 13 ) filename (3:6)='pl40' 

if ( iiu . eq . 12 ) filename (3:6)='pl20' 

if ( iiu . eq . 11 ) filename (3:6)='pl00’ 

if ( iiu . eq. 1 0 ) filename (3:6)='p090' 

if ( iiu . eq . 9) filename (3:6)='p080' 

if ( iiu . eq . 8 ) filename (3:6)='p070' 

if (iiu . eq. 7) filename (3:6)='p060' 

if ( iiu . eq . 6) filename (3:6)='p050' 

if ( iiu . eq . 5) filename (3:6)='p04Q' 

if (iiu .eq. 4) filename (3:6)='p030' 

if (iiu .eq. 3) filename ( 3 : 6) = ' p020 ' 

if (iiu .eq. 2 ) filename (3 : 6) = 'p010 ' 

if (iiu . eq. 1) filename (3:6)='p005' 

open (81, file=fi lename, status='old' ) 

icou=0 

do  100  i=l, 1000 

read (81, *, end=999) adum ( i) , tdum (i) 

100  icou=icou+l 

999  if (iiu.lt.ll)write (82, 110) iiu-1, icou, iiu-1, icou 

if (iiu.eq. 11) write (82,111) iiu-1, icou, iiu-1, icou 
if (iiu.gt.ll)write (82, 111) isep, icou, isep, icou 

110  format ( 6x,  'dimension  al ' , il ,  '  0  ( ' , i3,  ' ) , t ' , il ,  '  0 

111  formal (6x. 'dimension  al',i2, ' 0 ( ' , i3, ' ) , t ' , i2, '0 
if (iiu. It. 11) write (82, 115) iiu-1, icou 

if (iiu .eq. 11) write (82, 116) iiu-1, icou 
if (iiu.gt.ll)write (82, 116) isep, icou 

115  format ( 6x,  ' data  n  ' , il,  '0/ ' , i3,  '/ ' ) 

116  f ormat ( 6x, ' data  n',i2,'0/',i3,'/') 
do  200  i=l , icou 

if (iiu.lt . 1 1 ) write (82, 120) iiu-1, i, adum ( i ) , iiu-1 , 
if ( iiu . eq . 1 1 ) write (82, 121) iiu-1, i, adum ( i ) , iiu-1, 
if (iiu.gt.ll)write(82, 121) isep, i , adum ( i ) , isep, i, 


120 

format ( 6x, 

'data 

al ' 

,11, 

’ 0 ( ', i3,  ' )  / ' , el3 . 7, 

1  '/,t’,il 

,'0C, 

i3. 

'  )  /  ' 

, el3 . 7,  '/') 

121 

format ( 6x, 

'  data 

al ' 

,  12, 

' 0  (  ' , i 3 ,  ' ) / ' ,el3 .7, 

1  ' / , t ' , i2 

,  '0  (  ', 

13, 

'  )  /  ' 

, el 3 . 7 ,  '/') 

200 

cont inue 
close  (81) 

90 

continue 

stop 

end 


(  ’ , i3,  ' )  ’ ) 
(  ' , i3,  '  )  '  ) 


i , tdum ( i ) 
i , tdum ( i ) 
tdum  (  i ) 
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APPENDIX  D 

SAMPLE  OF  SUBROUTINE  TO  BE 
INCLUDED  INTO  A  ONE-,  TWO-,  OR  THREE-DIMENSIONAL 
HYDRODYNAMIC  COMPUTATION  TO  DETERMINE  A  NEW  BURN  FRACTION 


CDECK  USRDFB 

FUNCTION  USRDFB (FBRN, TIM, DT, PRS, SMP, NP, I , J, K, IMAX, JMAX, KMAX, 
1  SDETIN, PDETIN, PBREF, BUSR, CUSR, DUSR) 


C 

C 

C 

C 

C 

C 

C 

c 

c 


****************************************************************** 

***  THIS  IS  A  SAMPLE  USER-INPUT  FUNCTION  FOR  THE  EVALUATION 
OF  EXPLOSIVE  DECOMPOSITION 

***  THIS  MODULE  IS  ACCESSED  WHENEVER  BURN  MODEL  KBURN=-9  IS  CHOSEN 
FROM  INPUT 


c 

***  DEFINITION 

OF 

INPUT  PARAMETERS  - 

c 

USRDFB 

OUT 

CURRENT  VALUE  OF  BURN  RATE 

c 

FBRN 

INP 

'OLD'  VALUE  OF  BURN  FRACTION 

c 

OUT 

'NEW'  VALUE  OF  BURN  FRACTION 

c 

TIM 

INP 

CURRENT  TIME 

c 

DT 

INP 

CURRENT  TIME  STEP 

c 

PRS 

INP 

CURRENT  LOCAL  PRESSURE 

c 

SMP 

INP 

CURRENT  LOCAL  TEMPERATURE 

c 

NP 

INP 

CURRENT  EXTERNAL  MATERIAL  NUMBER 

c 

I 

INP 

CELL  INDEX  (I-DIRECTION) 

c 

J 

INP 

CELL  INDEX  ( J-DIRECTION) 

c 

K 

INP 

CELL  INDEX  (K-DIRECTION) 

c 

IMAX 

INP 

MAXIMUM  NUMBER  OF  CELLS 

(I-DIRECTION) 

c 

JMAX 

INP 

MAXIMUM  NUMBER  OF  CELLS 

(J-DIRECTION) 

c 

KMAX 

INP 

MAX T MUM  NUMBER  OF  CELLS 

(K-PIRECTION) 

c 

SDETIN 

INP 

USER  DEFINED  -  SUGGESTED 

INITIATION  TIME 

c 

PDETIN 

INP 

USER  DEFINED  -  SUGGESTED 

INITIATION  PRESSURE 

c 

PBREF 

INP 

USER  DEFINED  -  SUGGESTED 

REFERENCE  PRESSURE 

c 

BUSR 

INP 

USER  DEFINED  -  SUGGESTED 

NONE 

c 

CUSR 

INP 

USER  DEFINED  -  SUGGESTED 

NONE 

c 

DUSR 

INP 

USER  DEFINED  -  SUGGESTED 

NONE 

C 

C 

C 

C 

C 

C 

C 

C 

C 


THIS  ROUTINE  SHOULD  RETURN  - 

EITHER  A  BURN  RATE  BY  MEANS  OF  THE  FUNCTION  NAME 
OR  AN  UPDATED  VALUE  OF  THE  BURN  FRACTION  VIA  FBRN 

THE  APPROPRIATE  DEFINITION  IS  DONE  DEPENDING  ON  - 

IF  ON  OUTPUT  FBRN  IS  FOUND  TO  HAVE  CHANGED  ITS  VALUE, 
THEN  THE  UPDATED  VALUE  OF  FBRN  IS  USED 
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C  (IN  THIS  CASE  IT  IS  ADVISABLE  TO  DEFINE  USRDFB=0 . ) 

C  OTHERWISE  A  BURN  RATE  DEFINITION  IS  ASSUMED 

C 

Q  ★★★★★★★★***★★*★*******★**★★★****★****★★****★★*****★★*★*★*★*★**★*** 

IF (PRS . LT . 5 .E9) RETURN 

SFB=FBRN 

ZZFD=0 . 0 

IF (FBRN . GE . 1 . E-5 ) ZZFD=FBRN 

CALL  BEEN (ZZFD, TIM, DT, PRS, SMP,NP, I, J, K, IMAX, JMAX, KMAX, 

1  SDETIN, PDETIN, PBREF , BUSR, CUSR, DUSR) 

IF ( ZZFD. GE. 1 .E-5) THEN 
FBRN=ZZFD 
URSDFB=0 . 0 
ELSE 

FBRN=SFB 
URSDFB=1 . E-2  0 
ENDIF 
RETURN 
END 

SUBROUTINE  BEEN (FBRN, TIM, DT, PRS, SMP , NP , I , J, K, IMAX, JMAX, KMAX, 

1  SDETIN, PDETIN, PBREF, BUSR, CUSR, DUSR) 

DATA  CV/2.09E-5/,  AJWL/8 . 524011/,  BJWL/0 . 1801812/,  OMJWL/O.38/ 

DATA  R1JWL/4 . 6/,  R2JWL/1.3/,  CJWL/0 . 012066373/ 

DATA  TDT/ . 6666667/, TLIM/1 .E-3/ 

INCLUDE  ' LAMFILE . FOR ' 

INCLUDE  ' HSDATA . FOR  1 

common/ ic jckc/  IC (1000), JC (1000), KC (1000) ,FCC(1000) 

C  PUT  DT  INTO  MICROSECOND  UNITS  AND  SAVE  ORIGINAL  VALUE  OF  DT 
SAVEDT=DT 
DT=DT* 1 . E6 

C  IN  GENERAL,  PCJ  SHOULD  EQUAL  PBREF,  THE  CJ  PRESSURE  IN  THE  INPUT 
C  IE  PCJ  =  PBREF/ 1.E12  TO  GET  TO  MEGABARS. 

PCJ=PBREF/1 .E12 

C  PRES  IS  IN  UNITS  OF  MEGABARS,  PRS  IN  UNITS  OF  MICROBARS 
PRES=PRS/1 .E12 

C  USE  5  KBAR  AS  MINIMUM  PRESSURE  FOR  REACTION  TO  OCCUR  BASED  ON  GENERAL 
C  PURPOSE  PROGRAM  RESULTS 

C  MINIMUM  PRESSURE  SHOULD  CORRESPOND  TO  VALUE  IN  PO(NRS)  FROM  HS . INFO 
C  THE  POINT  BELOW  WHICH  THE  GENERAL  PURPOSE  PROGRAM  DOES  NOT  GO 
IF (PRES .LT. 0 . 00501) THEN 
DT=SAVEDT 
RETURN 
ENDIF 
IGOTOSUR=0 

IF (PRES .GT . 0 .200) IGOTOSUR=l 
C  DETERMINE  ALAMHS  BY  LOOKING  UP  TABLE  OF  DATA 

C  IF  FBRN  GE  AHS  USE  SURFACE  BURN,  IF  FBRN  LT  AHS  USE  NSWC  MODEL 
IREAD=0 
IPICK=0 

IF (BUSR . LT . 1 . E10 . AND . FBRN . LT . TLIM . AND . FBRN  . GE . 0 . 0 ) THEN 
IBUSR=IFIX(BUSR+0. 00001) 

DO  95  KDC=1, IBUSR 

IF ( IC (KDC) . EQ . I . AND . JC ( KDC ) . EQ . J . AND . KC (KDC) . EQ . K) THEN 
IPICK=1 

FBRN=FBRN+FCC (KDC) 

ENDIF 

95  CONTINUE 
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IREAD=1 

ENDIF 

SAVEFBRN=FBRN 
IF (IGOTOSUR.EQ. 1) GOTO  338 
NRSM 1 =NRS - 1 
DO  400  ICP=1 , NRSM1 

IF {PRES . GE . PO ( ICP+1 ) .AND. PRES. LT.PO(ICP) ) THEN 
ICP1=ICD+1 

AHS= (ALAMHS (ICP1) -ALAMHS (ICP) ) * (PRES-PO ( ICP ) ) / 

1  (PO (ICP1) -PO (ICP) ) +ALAMHS (ICP) 

TAHS= (THS (ICP1) -THS (ICP) ) * (PRES-PO ( ICP ) ) / 

1  (PO (ICP1) -PO (ICP) ) +THS (ICP) 

GOTO  401 
ENDIF 

400  CONTINUE 

401  CONTINUE 

IF (FBRN . LT . AHS . AND . FBRN . GE . 0 . 0 ) THEN 

C  USE  NSWC  MODEL  RELATIONS  TO  FIND  FBRN 
IF (PRES . GE . 0 . 1 8 ) THEN 

CALL  AFORP (AL18 0 , T1 80 , AL2 00 , T200 , N1 8 0 , N2 0 0 , FBRN, 

1  DT, PRES, PO (21 ) ,PO(l) , PCJ, THS (21 ) ,  THS (1) , AHS, TAHS) 

ELSEIF (PRES . LT . 0 . 18 . AND . PRES . GE . 0 . 16) THEN 
CALL  AFORP ( AL1 60 , T1 60 , AL18 0 , T180 , N1 60 , N1 8 0 , FBRN, 

1  DT , PRES , PO (41) , PO (21) , PCJ, THS (41) , THS (21) , AHS,  TAHS) 
ELSEIF (PRES .LT. 0 . 1 6 . AND . PRES . GE . 0 . 14) THEN 
CALL  AFORP ( AL1 40 , T1 40 , AL1 60 , T1 60 , N1 4 0 , N1 60 , FBRN, 

1  DT, PRES , PO ( 61 ) , PO (41) , PCJ, THS (61) , THS (41) , AHS, TAHS) 
ELSEIF (PRES . LT . 0 . 14 . AND . PRES .GE .0.12) THEN 
CALL  AFORP (AL120, T120, AL140, T140, N120, N140, FBRN, 

1  DT, PRES, PO (81) , PO ( 61) , PCJ, THS (81) , THS (61) , AHS, TAHS) 
ELSEIF (PRES . LT . 0 . 12 . AND . PRES . GE . 0 . 1 ) THEN 
CALL  AFORP ( AL100 , T10 0 , AL12 0 , T120 , N1 0 0 , N12 0 , FBRN, 

1  DT, PRES, PO (101) , PO (81) , PCJ, THS ( 10 1 )  ,  THS (81)  ,  AHS, TAHS) 
ELSEIF (PRES . LT . 0 . 1 . AND . PRES . GE . 0 . 0  9 ) THEN 
CALL  AFORP ( AL90 , T90 , ALIO 0 , T10 0 , N90 , N10 0 , FBRN, 

1  DT, PRES, PO( 111) ,PO(101) , PCJ, THS (111) , THS (101) , AHS, TAHS) 
ELSEIF (PRES . LT . 0 . 0 9 . AND . PRES . GE . 0 . 0 8 ) THEN 
CALL  AFORP ( AL8 0 , T80 , AL90 , T90 , N8 0 , N90 , FBRN, 

1  DT,  PRES,  PO(  121)  ,PO(lll)  ,  PCJ,  THS  (121.)  ,THS  (111)  ,  AHS,  TAHS) 
ELSEIF (PRES .LT. 0 . 08 . AND . PRES . GE . 0 . 07) THEN 
CALL  AFORP ( AL7 0 , T70 , AL8 0 , T8 0 , N7 0 , N80 , FBRN, 

1  DT, PRES, PO (131) , PO (121) , PCJ, THS ( 131 ) , THS (121) , AHS, TAHS) 
ELSEIF (PRES . LT . 0 . 07 . AND . PRES .GE .0.06) THEN 
CALL  AFORP ( AL60 , T60 , AL7 0 , T7 0 , N60 , N7 0 , FBRN, 

1  DT, PRES, PO (141) , PO ( 131 ) , PCJ, THS (141) , THS (131) , AHS, TAHS) 

ELSEIF (PRES . LT . 0 . 06 . AND . PRES . GE . 0 . 05) THEN 
CALL  AFORP (ALSO, T50, AL60, T60,N50,N60, FBRN, 

1  DT, PRES, PO (151) , PO (141) , PCJ, THS (151) , THS (141) , AHS, TAHS) 
ELSEIF (PRES . LT . 0 . 05 . AND . PRES .GE . 0 . 04 ) THEN 
CALL  AFORP ( AL40 , T4 0 , AL50 , T50 , N8 0 , N90 , FBRN, 

1  DT, PRES, PO (161) , PO (151) , PCJ, THS (161) , THS (151) , AHS, TAHS) 
ELSEIF (PRES. LT. 0.0 4 . AND . PRES . GE . 0 . 0 3 ) THEN 
CALL  AFORP ( AL30 , T30 , AL4 0 , T4 0 , N30 , N40 , FBRN, 

1  DT, PRES, PO (171) , PO (161) , PCJ, THS (171) , THS (161) , AHS, TAHS) 
ELSEIF (PRES . LT . 0 . 0 3 . AND . PRES . GE .0.02) THEN 
CALL  AFORP ( AL20 , T2 0 , AL30 , T30 , N20 , N3 0 , FBRN , 

1  DT,  PRES,  PO  (181)  ,  PO  (171'  ,  PC J,  THS  ( 1 8 1 )  ,  THS  (171)  ,  AHS,  TAHS) 
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ELSEIF (PRES . LT . 0 . 02 . AND . PRES . GE . 0 . 0 1 ) THEN 
CALL  AFORP ( ALIO , T10 , AL20 , T20 , N10 , N20 , FBRN, 

1  DT, PRES, PO (191) , PO (181) , PCJ, THS (191) , THS(181) , AHS, TAHS) 
ELSEIF (PRES. LT. 0.01) THEN 

CALL  AFORP (AL00 , TOO , ALIO, T10 , NOO , N1 0 , FBRN, 

1  DT , PRES , PO (NRS ) ,PO(191) , PCJ, THS (NRS ) , THS (191) , AHS, TAHS) 

ENDIF 

ENDIF 

AFORPFBRN=FBRN 
338  CONTINUE 

C  IF (FBRN . GE . AHS . AND . FBRN . LT . 1 . 0 ) THEN 

IF (FBRN . LT . 1 . 0 ) THEN 
FBRN=SAVEFBRN 

C  USE  SURFACE  BURN  RELATIONS  TO  FIND  FBRN 
CALL  SURKIM (PRES, FBRN, DT, AHS) 

ENDIF 

SURKIMFBRN=FBRN 

BRATE1= (AFORPFBRN-SAVEFBRN) /DT 
BRATE2= (SURKIMFBRN-SAVEFBRN) /DT 

C  FBRN  IS  CURRENTLY  FROM  SURKIM,  CHANGE  IF  BRATE1  FBRN  FROM  AFORP  IS 
LESS 

IF (BRATE1 . LE . BRATE2 ) FBRN=AFORPFBRN 
DT=SAVEDT 

IF (FBRN . GE . TLIM . AND . IREAD . EQ . 0 ) GOTO  774 
IF (FBRN . EQ . 0 . 0 . AND . IREAD . EQ . 0 ) GOTO  774 
IF (FBRN . LT . TLIM . AND . IREAD . EQ . 0 . AND . FBRN . GT . 0 . 0 ) THEN 
IF (BUSR . LT . 1 . E10 ) BUSR=BUSR+1 . 

IF (BUSR . GT . 1 . E10 ) BUSR=1 . 

IBUSR=IFIX (BUSR+0 .00001) 

IC ( IBUSR) =1 
JC ( IBUSR) =J 
KC (IBUSR) =K 
FCC( IBUSR) =FBRN 

ELSEIF (FBRN . GE . TLIM . AND . IREAD . EQ . 1 . AND . IPICK  FQ . 1 ) TnEN 
IKF=0 

DO  96  KDC=1, IBUSR 

IF(IKF.EQ.l) IC (KDC-1) =IC (KDC) 

IF  ( IKF  .  EQ  .  1 )  JC  ( KDC-1 )  =  JC  ( KDC') 

IF (IKF .EQ. 1) KC (KDC-1) =KC (KDC) 

IF (IKF.EQ. 1) FCC (KDC-1) =FCC (KDC) 

IF ( IC (KDC) . EQ . I . AND . JC (KDC) . EQ . J . AND . KC (KDC) . EQ . K) IKF=1 

96  CONTINUE 
IBUSR=IBUSR- 1 
BUSR=BUSR-1 . 

ELSEIF (FBRN . LT . TLIM . AND . IREAD . EQ . 1 . AND . IPICK . EQ . 1 ) THEN 
DO  97  KDC=1, IBUSR 

IF ( IC (KDC) . EQ . I . AND . JC ( KDC) . EQ . J . AND . KC (KDC) . EQ . K) 

1  FCC (KDC) =FBRN 

97  CONTINUE 

ELSEIF (FBRN . LT . TLIM .AND . IREAD . EQ . 1 . AND . IPICK . EQ . 0 ) THEN 
BUSR=BUSR+1 

IBUSR=IFIX (BUSR+0 .00001) 

IC (IBUSR) =1 
JC (IBUSR) =J 
KC( IBUSR) =K 
FCC (IBUSR) =FBRN 
ENDIF 
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774  CONTINUE 
RETURN 
END 

SUBROUTINE  AFORP ( A1 , T1 , A2 , T2 , N1 , N2 , FBRN, DT , PRES , PRES 1 , PRES2 , PC J, 
1  THS1, THS2, AHS, TAHS) 

FBRN:  INPUT  -OLD  BURN  FRACTION 

OUTPUT  -NEW  BURN  FRACTION 

DIMENSION  A1 (Nl) ,T1 (Nl) , A2 (N2) ,T2 (N2) , AAVE (1000) , TAVE (1000) 
INCLUDE  'HSDATA.FOR' 

IF (Nl . NE . N2 ) STOP  123 
QUANX= (PCJ/PRES) ** . 15 
QUANY= (PCJ/PRES) **.05 
N1M1=N1-1 
DO  100  IR=1 , Nl 

AAVE (IR) = (A2 (IR) -A1 (IR) ) * (PRES-PRES1) / (PRES2-PRES 1 ) +A1 (IR) 

TAVE (IR)  =  (T2 (IR) -T1 (IR) ) * (PRES-PRES1) / (PRES2-PRES1) +T1 (IR) 

100  CONTINUE 

IF (FBRN.EQ. 0 . 0) THEN 
TBN=DT 
ELSE 

FBNSCL=ALOG10 ( (FBRN/AHS ) **QUANY) 

C  NEED  TO  FIND  OLD  TIME 
DO  110  1=1 , N1M1 

IF (AAVE (I) . LT . FBNSCL . AND . AAVE ( I + 1 ) . GE . FBNSCL) THEN 
TBOSCL= (FBNSCL-AAVE (I) ) * (TAVE ( 1+ 1 ) -TAVE ( I ) ) / 

1  (AAVE (1+1) -AAVE (I) ) +TAVE (I) 

GOTO  111 
END  IF 

110  CONTINUE 

111  CONTINUE 

TBO=TAHS* ( (10 . *  *TBOSCL) ** (l./QUANX) ) 

TBN=TBO+DT 

ENDIF 

IF (TBN . GT . TAHS ) THEN 
FBRN=AHS 
USEDT=TAHS-TBO 
SDT=DT-USEDT 
IF ( SDT . GT . 0 . 0 ) DT=SDT 
RETURN 
ELSE 

TBNSCL=ALOG10 ( (TBN/TAHS ) **QUANX) 

DO  120  1=1 , N1M1 

IF (TAVE (I) . LT.TBNSCL. AND. TAVE (1+1) . GE . TBNSCL) THEN 
FBNSCL= (TBNSCL-TAVE (I) ) * ( AAVE ( I +1 ) -AAVE ( I ) ) / 

1  (TAVE (1+1) -TAVE (I) ) +AAVE (I) 

GOTO  121 
ENDIF 

120  CONTINUE 

121  CONTINUE 
ENDIF 

FBRN=AHS* ( (10 . **FBNSCL) ** (1 . /QUANY) ) 

RETURN 

END 

SUBROUTINE  SURKIM (PRES , FBRN, DT, AHS ) 

DATA  CV/2 .09E-5/,  AJWL/8 . 524011/,  BJWL/0 . 1801812/,  OMJWL/0.38/ 
DATA  R1JWL/4.6/,  R2JWL/1.3/,  CJWL/0 . 012066373/ 

DATA  TDT/ . 6 66 6 6 67 / , ALP / 1 . 2 / , BET/-0 . 4 /, DEL/ 32 500 . /,EP/2 .5/ 
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DATA  ROORG/ . Ol/,RIORG/ . 0 0 3 9 / , ANVAL/ 0 . 9/ 

SBR=3 . * . 00627/ROORG 

DSLAM=SBR* (FBRN* *TDT * ALP *PRES *  * (BET+ANVAL)  + 

1  (l.-FBRN) *  *TDT*DEL*PRES*  * (EP+ANVAL)  ) 

FBRN=FBRN+DT*DSLAM 
FBRN=AMIN1 (FBRN, 1 . ) 

RETURN 

END 
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